COVID19 UK

Prepare data

Read and format data

Prevalence

df_uk_prev <- read_csv('UK_timeseries_prep_2005.csv')

Parsed with column specification:
cols(
  ut_area = col_character(),
  date = col_character(),
  cumcase = col_double(),
  poptotal = col_double(),
  rate = col_double()
)

df_uk_prev <- df_uk_prev %>% 
  select(ut_area, date, rate) %>% 
  rename(rate_day = rate) %>%
  mutate(date = as.Date(date, "%d%b%Y"))



df_uk_prev2 <- read_csv('UK_timeseries_prep_0506_laua_all.csv')

Parsed with column specification:
cols(
  .default = col_double(),
  laua_name = col_character(),
  laua = col_character(),
  date = col_character()
)
See spec(...) for full column specifications.

df_uk_prev2 <- df_uk_prev2 %>% 
  select(laua, date, males, popdens, manufacturing, tourism, health,
         academic, medage, conservative, airport_dist, medinc, open,
         extra, agree, neuro, sci, rate) %>%
  rename(pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro,
         rate_day = rate) %>%
  mutate(date = as.Date(date, "%d%b%Y"))

df_uk_prev2 %>% select(-date, -rate_day) %>% 
  distinct() %>% write_csv('df_uk_pers_laua.csv')
  

Personality

df_uk_pers <- read_csv('timeseries_uk_utla_march9_april_09.csv')

Parsed with column specification:
cols(
  ut_area = col_character(),
  time = col_double(),
  areaname = col_character(),
  open = col_double(),
  extra = col_double(),
  agree = col_double(),
  neuro = col_double(),
  sci = col_double(),
  frequ = col_double(),
  ut_name = col_character(),
  poptotal = col_double(),
  rate_day = col_double()
)

df_uk_pers <- df_uk_pers %>% 
  select(ut_area, open, agree, neuro, sci, extra) %>% 
  dplyr::rename(pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro) %>%
  distinct()

df_uk_pers %>% write_csv('df_uk_pers_ut.csv')

df_uk_pers_nuts <- read_csv('UK_socdist_fb_nuts3.csv')

Parsed with column specification:
cols(
  nuts3 = col_character(),
  date = col_date(format = ""),
  all_day_bing_tiles_visited_relat = col_double(),
  all_day_ratio_single_tile_users = col_double(),
  open = col_double(),
  extra = col_double(),
  agree = col_double(),
  neuro = col_double(),
  sci = col_double(),
  frequ = col_double(),
  nuts3_name = col_character(),
  runday = col_double()
)

df_uk_pers_nuts <- df_uk_pers_nuts %>% 
  select(nuts3, open, sci, extra, agree, neuro) %>%
  dplyr::rename(pers_o = open, 
                pers_c = sci,
                pers_e = extra,
                pers_a = agree,
                pers_n = neuro) %>%
  distinct()

df_uk_pers_nuts 

NA

Social distancing

fb_files <- list.files('../FB Data/UK individual files',
                       '*.csv', full.names = T)

df_uk_socdist <- fb_files %>% 
  map(read.csv) %>% bind_rows()

df_uk_socdist <- df_uk_socdist %>%
  select(ds, all_day_bing_tiles_visited_relative_change,
         all_day_ratio_single_tile_users, external_polygon_id) %>%
  rename(date = ds,
         nuts3 = external_polygon_id,
         socdist_tiles = all_day_bing_tiles_visited_relative_change,
         socdist_single_tile = all_day_ratio_single_tile_users) %>%
  mutate(nuts3 = as.character(nuts3),
         date = as.Date(date)) %>%
  arrange(nuts3, date) %>%
  drop_na()

df_uk_socdist %>% select(nuts3) %>% distinct() %>% nrow()

Controls

df_uk_ctrl_nuts <- read_csv("controls_UK_nuts3.csv")

Parsed with column specification:
cols(
  nuts3 = col_character(),
  nuts3_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_nuts <- df_uk_ctrl_nuts %>% select(-nuts3_name)
df_uk_ctrl_nuts

df_uk_ctrl_ut <- read_csv("controls_UK_ut.csv")

Parsed with column specification:
cols(
  ut_area = col_character(),
  ut_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_ut <- df_uk_ctrl_ut %>% select(-ut_name)
df_uk_ctrl_ut

NA
NA

Merge prevalence data

# create sequence of dates
date_sequence <- seq.Date(min(df_uk_prev2$date),
                          max(df_uk_prev2$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge day index with gps data
df_uk_prev2 = df_uk_prev2 %>% 
  merge(df_dates, by='date') %>% 
  arrange(laua) %>%
  as_tibble()

df_uk_prev2

NA

Merge social distancing data

Check timeframes

df_uk_prev$date %>% summary()
df_uk_socdist$date %>% summary()

Control for weekend effect in social distancing

easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)


df_uk_loess <- df_uk_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$nuts3) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_uk_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'nuts3', value = 'loess') %>% 
  group_by(nuts3) %>% 
  mutate(time = row_number())

df_uk_loess_2 <- df_uk_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$nuts3) %>%
  map(~ loess(socdist_tiles ~ time, data = .)) %>%
  map(predict, 1:max(df_uk_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'nuts3', value = 'loess') %>%
  rename(loess_2 = loess) %>%
  group_by(nuts3) %>% 
  mutate(time = row_number())

df_uk_socdist <- df_uk_socdist %>% 
  merge(df_uk_loess, by=c('nuts3', 'time')) %>% 
  merge(df_uk_loess_2, by=c('nuts3', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess, socdist_single_tile),
         socdist_tiles_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess_2, socdist_tiles)) %>%
  arrange(nuts3, time) %>% 
  select(-weekday)

df_uk_socdist <- df_uk_socdist %>% drop_na() %>% mutate(time = time-1)

Explore data

Plot prevalence over time

df_uk_prev %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_prev %>% mutate(prev_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_uk_prev2 %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=laua, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_prev2 %>% mutate(prev_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=laua, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Plot social distancing over time

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(socdist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~socdist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_tiles)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(socdist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_tiles)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~socdist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_uk_socdist <- df_uk_socdist %>% 
  mutate(socdist_single_tile = socdist_single_tile_clean,
         socdist_tiles = socdist_tiles_clean) %>% 
  select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)

Correlations

df_uk_prev %>% group_by(ut_area) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-ut_area, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3) %>%
  as.data.frame()

df_uk_socdist %>% group_by(nuts3) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-nuts3, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3) %>% 
  as.data.frame()

NA

Rescale Data

lvl2_scaled_ut <- df_uk_prev %>% 
  dplyr::select(-time, -date, -rate_day) %>% 
  distinct() %>% 
  mutate_at(vars(-ut_area), scale)

lvl1_scaled_ut <- df_uk_prev %>% select(ut_area, time, rate_day)

df_uk_prev_scaled <- plyr::join(lvl1_scaled_ut, lvl2_scaled_ut, by = 'ut_area')

df_uk_prev_scaled

lvl2_scaled_laua <- df_uk_prev2 %>% 
  dplyr::select(-time, -date, -rate_day) %>% 
  distinct() %>% 
  mutate_at(vars(-laua), scale)

lvl1_scaled_laua <- df_uk_prev2 %>% select(laua, time, rate_day)

df_uk_prev2_scaled <- plyr::join(lvl1_scaled_laua, lvl2_scaled_laua, by = 'laua')

df_uk_prev2_scaled

NA

lvl2_scaled_nuts <- df_uk_socdist %>% 
  dplyr::select(-time, -date, -socdist_tiles, -socdist_single_tile) %>% 
  distinct() %>% 
  mutate_at(vars(-nuts3), scale)

lvl1_scaled_nuts <- df_uk_socdist %>% 
  select(nuts3, time, socdist_single_tile, socdist_tiles) %>% 
  mutate_at(vars(-nuts3, -time), scale)

df_uk_socdist_scaled <- plyr::join(lvl1_scaled_nuts, lvl2_scaled_nuts, by = 'nuts3')

df_uk_socdist_scaled

NA

Predict Prevalence UT Level

Extract first day of covid outbreak

# get onset day
df_uk_onset_prev <- df_uk_prev_scaled %>% 
  group_by(ut_area) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time))
  
# merge with county data
df_uk_onset_prev <- df_uk_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  left_join(df_uk_onset_prev, by = 'ut_area')

# handle censored data
df_uk_onset_prev <- df_uk_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_uk_prev$date)))+1))

Extract slopes

# cut time series before onset
df_uk_prev_scaled <- df_uk_prev_scaled %>% 
  group_by(ut_area) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_uk_prev_scaled <- df_uk_prev_scaled %>%
  group_by(ut_area) %>%
  filter(n() == 30) %>%
  ungroup()

# log transform prevalence data 
df_uk_prev_scaled <- df_uk_prev_scaled %>% 
  mutate(rate_day = log(rate_day))

# extract slope prevalence
df_uk_slope_prev <- df_uk_prev_scaled %>% split(.$ut_area) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('ut_area') %>% 
  rename(slope_prev = '.')

# merge with county data
df_uk_slope_prev <- df_uk_onset_prev %>% 
  inner_join(df_uk_slope_prev, by = 'ut_area') %>%
  drop_na()

# standardize slopes
df_uk_slope_prev <- df_uk_slope_prev %>% 
  mutate(slope_prev = scale(slope_prev),
         onset_prev = scale(onset_prev))

Explore distributions

df_uk_onset_prev %>% ggplot(aes(onset_prev)) + geom_histogram()
df_uk_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram()

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_uk_onset_prev)
cox_onset_prev %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_uk_onset_prev)

  n= 149, number of events= 149 

            coef exp(coef)  se(coef)      z Pr(>|z|)  
pers_o  0.322264  1.380249  0.166864  1.931   0.0534 .
pers_c  0.036475  1.037149  0.135119  0.270   0.7872  
pers_e -0.013895  0.986201  0.143396 -0.097   0.9228  
pers_a  0.005323  1.005338  0.106927  0.050   0.9603  
pers_n -0.092400  0.911740  0.106197 -0.870   0.3843  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.3802     0.7245    0.9952     1.914
pers_c    1.0371     0.9642    0.7958     1.352
pers_e    0.9862     1.0140    0.7446     1.306
pers_a    1.0053     0.9947    0.8153     1.240
pers_n    0.9117     1.0968    0.7404     1.123

Concordance= 0.605  (se = 0.027 )
Likelihood ratio test= 12.64  on 5 df,   p=0.03
Wald test            = 13.47  on 5 df,   p=0.02
Score (logrank) test = 13.78  on 5 df,   p=0.02

# predict onset from personality with controls
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative,
                             data = df_uk_onset_prev)
cox_onset_prev_ctrl %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative, 
    data = df_uk_onset_prev)

  n= 144, number of events= 144 
   (5 observations deleted due to missingness)

                  coef exp(coef) se(coef)      z Pr(>|z|)  
pers_o         0.13431   1.14375  0.26537  0.506   0.6128  
pers_c         0.27571   1.31746  0.21410  1.288   0.1978  
pers_e        -0.28546   0.75167  0.16166 -1.766   0.0774 .
pers_a         0.15618   1.16904  0.12936  1.207   0.2273  
pers_n        -0.04289   0.95801  0.13821 -0.310   0.7563  
airport_dist   0.08819   1.09219  0.10241  0.861   0.3892  
males         -0.28751   0.75013  0.11525 -2.495   0.0126 *
popdens        0.31764   1.37388  0.22294  1.425   0.1542  
manufacturing -0.01434   0.98576  0.15053 -0.095   0.9241  
tourism        0.20101   1.22263  0.12897  1.559   0.1191  
health        -0.13046   0.87769  0.09930 -1.314   0.1889  
academic       0.41323   1.51170  0.25879  1.597   0.1103  
medinc        -0.10312   0.90201  0.15281 -0.675   0.4998  
medage        -0.54104   0.58214  0.21413 -2.527   0.0115 *
conservative  -0.05175   0.94957  0.28447 -0.182   0.8557  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.1438     0.8743    0.6799    1.9241
pers_c           1.3175     0.7590    0.8660    2.0044
pers_e           0.7517     1.3304    0.5475    1.0319
pers_a           1.1690     0.8554    0.9072    1.5064
pers_n           0.9580     1.0438    0.7307    1.2561
airport_dist     1.0922     0.9156    0.8936    1.3350
males            0.7501     1.3331    0.5985    0.9402
popdens          1.3739     0.7279    0.8875    2.1267
manufacturing    0.9858     1.0144    0.7339    1.3240
tourism          1.2226     0.8179    0.9496    1.5742
health           0.8777     1.1394    0.7225    1.0663
academic         1.5117     0.6615    0.9103    2.5104
medinc           0.9020     1.1086    0.6686    1.2170
medage           0.5821     1.7178    0.3826    0.8857
conservative     0.9496     1.0531    0.5437    1.6583

Concordance= 0.643  (se = 0.027 )
Likelihood ratio test= 37.26  on 15 df,   p=0.001
Wald test            = 39.13  on 15 df,   p=6e-04
Score (logrank) test = 40.98  on 15 df,   p=3e-04

Predict prevalence slopes with linear models

# predict slopes from personality
lm_slope_prev <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_uk_slope_prev)
lm_slope_prev %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_uk_slope_prev)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7474 -0.3538  0.2157  0.5971  1.6747 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.003326   0.081259   0.041 0.967406    
pers_o      -0.335161   0.160109  -2.093 0.038150 *  
pers_c      -0.471364   0.134224  -3.512 0.000602 ***
pers_e       0.045912   0.149655   0.307 0.759467    
pers_a       0.087439   0.115589   0.756 0.450657    
pers_n      -0.172004   0.120288  -1.430 0.154997    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9747 on 138 degrees of freedom
Multiple R-squared:  0.08317,   Adjusted R-squared:  0.04995 
F-statistic: 2.504 on 5 and 138 DF,  p-value: 0.03328

lm_slope_prev %>% confint(level=0.9)

                   5 %        95 %
(Intercept) -0.1312354  0.13788821
pers_o      -0.6002971 -0.07002450
pers_c      -0.6936350 -0.24909386
pers_e      -0.2019112  0.29373593
pers_a      -0.1039721  0.27885073
pers_n      -0.3711981  0.02719003

# predict slopes from personality with controls
lm_slope_prev_ctrl <- lm(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev,
                         data = df_uk_slope_prev)
lm_slope_prev_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative + 
    onset_prev, data = df_uk_slope_prev)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.05366 -0.46070  0.00142  0.46146  1.97316 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.006093   0.062803   0.097   0.9229    
pers_o        -0.107115   0.189739  -0.565   0.5734    
pers_c        -0.314027   0.148365  -2.117   0.0362 *  
pers_e        -0.132175   0.124631  -1.061   0.2909    
pers_a         0.244855   0.109204   2.242   0.0267 *  
pers_n        -0.172078   0.111285  -1.546   0.1245    
airport_dist  -0.127434   0.082936  -1.537   0.1269    
males         -0.220415   0.094095  -2.342   0.0207 *  
popdens        0.313726   0.160336   1.957   0.0526 .  
manufacturing  0.077796   0.107872   0.721   0.4721    
tourism       -0.016370   0.091325  -0.179   0.8580    
health        -0.033028   0.079300  -0.416   0.6778    
academic       0.125265   0.192659   0.650   0.5167    
medinc         0.146443   0.116369   1.258   0.2105    
medage        -0.059607   0.155676  -0.383   0.7024    
conservative   0.134361   0.208331   0.645   0.5201    
onset_prev     0.589208   0.069384   8.492  4.6e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7522 on 127 degrees of freedom
Multiple R-squared:  0.4975,    Adjusted R-squared:  0.4342 
F-statistic: 7.859 on 16 and 127 DF,  p-value: 1.121e-12

lm_slope_prev_ctrl %>% confint(level=0.9)

                      5 %         95 %
(Intercept)   -0.09796849  0.110153915
pers_o        -0.42150130  0.207271247
pers_c        -0.55985809 -0.068195391
pers_e        -0.33868197  0.074331811
pers_a         0.06391090  0.425799311
pers_n        -0.35647153  0.012315179
airport_dist  -0.26485443  0.009986441
males         -0.37632416 -0.064506045
popdens        0.04805812  0.579393990
manufacturing -0.10094080  0.256532733
tourism       -0.16768973  0.134950274
health        -0.16442352  0.098367638
academic      -0.19395953  0.444488789
medinc        -0.04637384  0.339260214
medage        -0.31755259  0.198338423
conservative  -0.21083119  0.479552475
onset_prev     0.47424351  0.704172203

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev,
                           data = df_uk_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp

       pers_o        pers_c        pers_e        pers_a        pers_n  airport_dist 
 0.0128740502  0.0185184287 -0.0018344699 -0.0034459553 -0.0047870867  0.0471837405 
        males       popdens manufacturing       tourism        health      academic 
-0.0009166694  0.0159040281 -0.0020926242  0.0182207566 -0.0015815550  0.0023041324 
       medinc        medage  conservative    onset_prev 
 0.0124853594  0.0028673495  0.0002912342  0.5122363597 

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond

       pers_o        pers_c        pers_e        pers_a        pers_n  airport_dist 
  0.013531981   0.015370914   0.001291142  -0.004426900  -0.002525750   0.053357454 
        males       popdens manufacturing       tourism        health      academic 
  0.002060008   0.018857165  -0.003769689   0.013425660   0.000157888   0.004189666 
       medinc        medage  conservative    onset_prev 
  0.013498891   0.007233071   0.001083141   0.443129295 

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Prevalence LAD Level

Extract first day of covid outbreak

# get onset day
df_uk_onset_prev_laua <- df_uk_prev2_scaled %>% 
  group_by(laua) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time))
  
# merge with county data
df_uk_onset_prev_laua <- df_uk_prev2_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  left_join(df_uk_onset_prev_laua, by = 'laua')

# handle censored data
df_uk_onset_prev_laua <- df_uk_onset_prev_laua %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_uk_prev2$date)))+1))

Extract slopes

# cut time series before onset
df_uk_prev2_scaled <- df_uk_prev2_scaled %>% 
  group_by(laua) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_uk_prev2_scaled <- df_uk_prev2_scaled %>%
  group_by(laua) %>%
  filter(n() == 30) %>%
  ungroup()

# log transform prevalence data 
df_uk_prev2_scaled <- df_uk_prev2_scaled %>% 
  mutate(rate_day = log(rate_day))

# extract slope prevalence
df_uk_slope_prev_laua <- df_uk_prev2_scaled %>% split(.$laua) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('laua') %>% 
  rename(slope_prev = '.')

# merge with county data
df_uk_slope_prev_laua <- df_uk_onset_prev_laua %>% 
  inner_join(df_uk_slope_prev_laua, by = 'laua') %>%
  drop_na()

# standardize slopes
df_uk_slope_prev_laua <- df_uk_slope_prev_laua %>% 
  mutate(slope_prev = scale(slope_prev),
         onset_prev = scale(onset_prev))

Explore distributions

df_uk_onset_prev_laua %>% ggplot(aes(onset_prev)) + geom_histogram()
df_uk_slope_prev_laua %>% ggplot(aes(slope_prev)) + geom_histogram()

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev_laua <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_uk_onset_prev_laua)
cox_onset_prev_laua %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_uk_onset_prev_laua)

  n= 312, number of events= 312 

           coef exp(coef) se(coef)      z Pr(>|z|)   
pers_o  0.11799   1.12524  0.07472  1.579   0.1143   
pers_c -0.17400   0.84030  0.07851 -2.216   0.0267 * 
pers_e  0.20837   1.23167  0.07505  2.776   0.0055 **
pers_a  0.05375   1.05522  0.06484  0.829   0.4071   
pers_n  0.01251   1.01259  0.07158  0.175   0.8613   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.1252     0.8887    0.9720    1.3027
pers_c    0.8403     1.1901    0.7205    0.9801
pers_e    1.2317     0.8119    1.0632    1.4269
pers_a    1.0552     0.9477    0.9293    1.1982
pers_n    1.0126     0.9876    0.8800    1.1651

Concordance= 0.602  (se = 0.02 )
Likelihood ratio test= 27.18  on 5 df,   p=5e-05
Wald test            = 29.19  on 5 df,   p=2e-05
Score (logrank) test = 29.1  on 5 df,   p=2e-05

# predict onset from personality with controls
cox_onset_prev_laua_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative,
                             data = df_uk_onset_prev_laua)
cox_onset_prev_laua_ctrl %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative, 
    data = df_uk_onset_prev_laua)

  n= 312, number of events= 312 

                    coef  exp(coef)   se(coef)      z Pr(>|z|)    
pers_o        -0.0341890  0.9663888  0.1137888 -0.300 0.763826    
pers_c        -0.0376070  0.9630914  0.1059875 -0.355 0.722721    
pers_e         0.0811992  1.0845869  0.0807194  1.006 0.314443    
pers_a         0.0133296  1.0134188  0.0815961  0.163 0.870234    
pers_n        -0.0004303  0.9995697  0.0781902 -0.006 0.995609    
airport_dist  -0.0295206  0.9709108  0.0739004 -0.399 0.689550    
males         -0.2387449  0.7876157  0.0701625 -3.403 0.000667 ***
popdens        0.2236609  1.2506468  0.1028405  2.175 0.029643 *  
manufacturing -0.1338256  0.8747426  0.0770596 -1.737 0.082449 .  
tourism        0.0758201  1.0787685  0.0752577  1.007 0.313707    
health        -0.0853930  0.9181514  0.0624438 -1.368 0.171463    
academic       0.1458217  1.1569899  0.1170275  1.246 0.212747    
medinc        -0.1750250  0.8394361  0.0913846 -1.915 0.055460 .  
medage        -0.2801971  0.7556348  0.1237521 -2.264 0.023563 *  
conservative  -0.1640329  0.8487141  0.1402000 -1.170 0.242004    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           0.9664     1.0348    0.7732    1.2078
pers_c           0.9631     1.0383    0.7824    1.1855
pers_e           1.0846     0.9220    0.9259    1.2705
pers_a           1.0134     0.9868    0.8636    1.1892
pers_n           0.9996     1.0004    0.8575    1.1651
airport_dist     0.9709     1.0300    0.8400    1.1222
males            0.7876     1.2697    0.6864    0.9037
popdens          1.2506     0.7996    1.0223    1.5299
manufacturing    0.8747     1.1432    0.7521    1.0174
tourism          1.0788     0.9270    0.9308    1.2502
health           0.9182     1.0891    0.8124    1.0377
academic         1.1570     0.8643    0.9198    1.4553
medinc           0.8394     1.1913    0.7018    1.0041
medage           0.7556     1.3234    0.5929    0.9631
conservative     0.8487     1.1783    0.6448    1.1171

Concordance= 0.649  (se = 0.018 )
Likelihood ratio test= 71.42  on 15 df,   p=2e-09
Wald test            = 81.62  on 15 df,   p=4e-11
Score (logrank) test = 87.17  on 15 df,   p=3e-12

Predict prevalence slopes with linear models

# predict slopes from personality
lm_slope_prev_laua <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_uk_slope_prev_laua)
lm_slope_prev_laua %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_uk_slope_prev_laua)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2133 -0.4129  0.1032  0.5263  2.2010 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.413e-16  5.283e-02   0.000  1.00000    
pers_o      -2.208e-01  7.310e-02  -3.020  0.00274 ** 
pers_c      -4.870e-01  7.765e-02  -6.272 1.22e-09 ***
pers_e       8.785e-02  7.165e-02   1.226  0.22108    
pers_a       4.990e-02  7.034e-02   0.709  0.47859    
pers_n      -7.062e-02  7.215e-02  -0.979  0.32846    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9331 on 306 degrees of freedom
Multiple R-squared:  0.1432,    Adjusted R-squared:  0.1292 
F-statistic: 10.23 on 5 and 306 DF,  p-value: 4.425e-09

lm_slope_prev_laua %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.08716012  0.08716012
pers_o      -0.34138810 -0.10016782
pers_c      -0.61509638 -0.35888294
pers_e      -0.03035627  0.20605391
pers_a      -0.06615124  0.16595921
pers_n      -0.18966710  0.04841944

# predict slopes from personality with controls
lm_slope_prev_laua_ctrl <- lm(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev,
                         data = df_uk_slope_prev_laua)
lm_slope_prev_laua_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative + 
    onset_prev, data = df_uk_slope_prev_laua)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.36587 -0.47934  0.07511  0.48233  2.05828 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.157e-15  4.388e-02   0.000 1.000000    
pers_o        -1.082e-01  8.889e-02  -1.217 0.224581    
pers_c        -3.216e-01  8.221e-02  -3.912 0.000114 ***
pers_e         5.252e-02  6.319e-02   0.831 0.406543    
pers_a         5.110e-02  6.473e-02   0.789 0.430500    
pers_n        -6.529e-02  6.455e-02  -1.012 0.312594    
airport_dist  -1.143e-01  5.733e-02  -1.994 0.047023 *  
males         -2.664e-01  5.906e-02  -4.510 9.36e-06 ***
popdens        1.996e-01  8.221e-02   2.428 0.015780 *  
manufacturing -6.376e-02  5.946e-02  -1.072 0.284432    
tourism       -6.622e-04  6.027e-02  -0.011 0.991240    
health         3.748e-02  5.184e-02   0.723 0.470214    
academic       5.999e-02  9.881e-02   0.607 0.544258    
medinc        -8.994e-02  6.897e-02  -1.304 0.193197    
medage        -2.384e-01  9.994e-02  -2.386 0.017679 *  
conservative   7.081e-02  1.133e-01   0.625 0.532461    
onset_prev     5.111e-01  4.750e-02  10.761  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7751 on 295 degrees of freedom
Multiple R-squared:  0.4301,    Adjusted R-squared:  0.3992 
F-statistic: 13.91 on 16 and 295 DF,  p-value: < 2.2e-16

lm_slope_prev_laua_ctrl %>% confint(level=0.9)

                      5 %        95 %
(Intercept)   -0.07240974  0.07240974
pers_o        -0.25485049  0.03849370
pers_c        -0.45728147 -0.18596760
pers_e        -0.05173975  0.15677729
pers_a        -0.05570562  0.15790033
pers_n        -0.17179645  0.04121405
airport_dist  -0.20895058 -0.01974658
males         -0.36381053 -0.16891665
popdens        0.06395855  0.33526073
manufacturing -0.16187386  0.03434819
tourism       -0.10010886  0.09878437
health        -0.04805089  0.12301281
academic      -0.10305098  0.22302120
medinc        -0.20374181  0.02385425
medage        -0.40333985 -0.07352087
conservative  -0.11614187  0.25776935
onset_prev     0.43277274  0.58952577

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev_laua <- cforest(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev,
                           data = df_uk_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev_laua, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev_laua, conditional = T, nperm = 1)

crf_slope_prev_varimp

       pers_o        pers_c        pers_e        pers_a        pers_n  airport_dist 
 6.004034e-03  2.052057e-02 -2.477124e-03 -5.399878e-03 -1.560624e-03  5.323515e-02 
        males       popdens manufacturing       tourism        health      academic 
-1.243299e-03  1.419716e-02  3.355199e-04  1.913226e-02 -5.297063e-04 -1.180278e-03 
       medinc        medage  conservative    onset_prev 
 1.182686e-02  2.453000e-03 -1.151113e-05  5.313719e-01 

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond

       pers_o        pers_c        pers_e        pers_a        pers_n  airport_dist 
 0.0081269597  0.0162285119  0.0002770556 -0.0032828696 -0.0007791452  0.0492546778 
        males       popdens manufacturing       tourism        health      academic 
-0.0012493988  0.0162960501  0.0026788967  0.0173285064  0.0030926219  0.0048356668 
       medinc        medage  conservative    onset_prev 
 0.0095891577  0.0054089674  0.0018216539  0.4410251062 

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Social Distancing

Change point analysis

# keep only counties with full data
nuts_complete <- df_uk_socdist_scaled %>% 
  group_by(nuts3) %>% 
  summarize(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$nuts3

# run changepoint analysis
df_uk_socdist_cpt_results <- df_uk_socdist_scaled %>% 
  select(nuts3, socdist_single_tile) %>%
  filter(nuts3 %in% nuts_complete) %>% 
  split(.$nuts3) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

df_uk_socdist_cpt_results_2 <- df_uk_socdist_scaled %>% 
  select(nuts3, socdist_tiles) %>%
  filter(nuts3 %in% nuts_complete) %>% 
  split(.$nuts3) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_tiles),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_uk_socdist_cpt_day <- df_uk_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('nuts3')

df_uk_socdist_cpt_day_2 <- df_uk_socdist_cpt_results_2 %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist_2 = '.') %>%
  rownames_to_column('nuts3')

# calculate mean differences
df_uk_socdist_cpt_mean_diff <- df_uk_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('nuts3')

df_uk_socdist_cpt_mean_diff_2 <- df_uk_socdist_cpt_results_2 %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ -.[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist_2 = '.') %>%
  rownames_to_column('nuts3')

# calculate means 
df_uk_socdist_mean <- df_uk_socdist_scaled %>% 
  group_by(nuts3) %>%
  summarize(mean_socdist = mean(socdist_single_tile))

df_uk_socdist_mean_2 <- df_uk_socdist_scaled %>% 
  group_by(nuts3) %>%
  summarize(mean_socdist_2 = -mean(socdist_tiles))

# merge with county data
nuts_ut_key <- read_csv('nuts3_ut.csv')

Parsed with column specification:
cols(
  nuts3 = col_character(),
  ut_area = col_character()
)

df_uk_cpt_socdist <- df_uk_socdist_scaled %>% 
  select(-time, -socdist_single_tile, -socdist_tiles) %>%
  distinct() %>% 
  left_join(df_uk_socdist_cpt_day, by='nuts3') %>%
  left_join(df_uk_socdist_cpt_day_2, by='nuts3') %>%
  left_join(df_uk_socdist_cpt_mean_diff, by='nuts3') %>%
  left_join(df_uk_socdist_cpt_mean_diff_2, by='nuts3') %>%
  left_join(df_ger_socdist_mean, by='nuts3') %>%
  left_join(df_ger_socdist_mean_2, by='nuts3') %>%
  left_join(nuts_ut_key, by='nuts3') %>% 
  left_join(select(df_uk_onset_prev, ut_area, onset_prev), by='ut_area') %>%
  left_join(select(df_uk_slope_prev, ut_area, slope_prev), by='ut_area') %>%
  select(-ut_area)

Error in tbl_vars_dispatch(x) : object 'df_ger_socdist_mean' not found

df_uk_cpt_socdist$cpt_day_socdist %>% hist()

df_uk_cpt_socdist$mean_diff_socdist %>% hist()

df_uk_cpt_socdist$mean_socdist %>% hist()

df_uk_cpt_socdist$cpt_day_socdist_2 %>% hist()

df_uk_cpt_socdist$mean_diff_socdist_2 %>% hist()

df_uk_cpt_socdist$mean_socdist_2 %>% hist()

cor(df_uk_cpt_socdist$mean_diff_socdist, 
    df_uk_cpt_socdist$mean_socdist)

          [,1]
[1,] 0.9756932

cor(df_uk_cpt_socdist$mean_diff_socdist_2, 
    df_uk_cpt_socdist$mean_socdist_2)

          [,1]
[1,] 0.9855565

for(i in head(df_uk_socdist_cpt_results_2, 5)){
  plot(i)
}

NA

for(i in head(df_uk_socdist_cpt_results_2, 5)){
  plot(i)
}

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist <- coxph(Surv(cpt_day_socdist, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_uk_cpt_socdist)
cox_cpt_socdist %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_uk_cpt_socdist)

  n= 131, number of events= 131 

            coef exp(coef)  se(coef)      z Pr(>|z|)
pers_o  0.037677  1.038396  0.159307  0.237    0.813
pers_c  0.006832  1.006855  0.144736  0.047    0.962
pers_e  0.002145  1.002148  0.160157  0.013    0.989
pers_a -0.054244  0.947201  0.134632 -0.403    0.687
pers_n -0.034435  0.966152  0.130787 -0.263    0.792

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.0384     0.9630    0.7599     1.419
pers_c    1.0069     0.9932    0.7582     1.337
pers_e    1.0021     0.9979    0.7322     1.372
pers_a    0.9472     1.0557    0.7275     1.233
pers_n    0.9662     1.0350    0.7477     1.248

Concordance= 0.941  (se = 0.033 )
Likelihood ratio test= 0.82  on 5 df,   p=1
Wald test            = 0.86  on 5 df,   p=1
Score (logrank) test = 0.86  on 5 df,   p=1

# predict hazard from personality with controls
cox_cpt_socdist_ctrl <- coxph(Surv(cpt_day_socdist, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                  data = df_uk_cpt_socdist)
cox_cpt_socdist_ctrl %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_dist + males + popdens + 
    manufacturing + tourism + health + academic + medinc + medage + 
    conservative + onset_prev + slope_prev, data = df_uk_cpt_socdist)

  n= 130, number of events= 130 
   (1 observation deleted due to missingness)

                   coef exp(coef)  se(coef)      z Pr(>|z|)
pers_o        -0.054009  0.947424  0.259655 -0.208    0.835
pers_c         0.015151  1.015266  0.207295  0.073    0.942
pers_e        -0.011907  0.988164  0.169616 -0.070    0.944
pers_a        -0.038884  0.961862  0.168035 -0.231    0.817
pers_n        -0.036467  0.964189  0.145252 -0.251    0.802
airport_dist   0.027010  1.027378  0.124971  0.216    0.829
males         -0.065868  0.936255  0.163367 -0.403    0.687
popdens        0.133712  1.143063  0.239108  0.559    0.576
manufacturing  0.052628  1.054038  0.140610  0.374    0.708
tourism        0.017427  1.017579  0.131044  0.133    0.894
health        -0.026273  0.974069  0.119258 -0.220    0.826
academic       0.055605  1.057180  0.267586  0.208    0.835
medinc        -0.049497  0.951708  0.168257 -0.294    0.769
medage        -0.051444  0.949857  0.232922 -0.221    0.825
conservative  -0.052351  0.948996  0.288138 -0.182    0.856
onset_prev     0.001703  1.001704  0.015047  0.113    0.910
slope_prev    -0.021361  0.978865  0.126324 -0.169    0.866

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           0.9474     1.0555    0.5695     1.576
pers_c           1.0153     0.9850    0.6763     1.524
pers_e           0.9882     1.0120    0.7087     1.378
pers_a           0.9619     1.0397    0.6920     1.337
pers_n           0.9642     1.0371    0.7253     1.282
airport_dist     1.0274     0.9734    0.8042     1.313
males            0.9363     1.0681    0.6797     1.290
popdens          1.1431     0.8748    0.7154     1.826
manufacturing    1.0540     0.9487    0.8001     1.388
tourism          1.0176     0.9827    0.7871     1.316
health           0.9741     1.0266    0.7710     1.231
academic         1.0572     0.9459    0.6257     1.786
medinc           0.9517     1.0507    0.6844     1.324
medage           0.9499     1.0528    0.6017     1.499
conservative     0.9490     1.0537    0.5395     1.669
onset_prev       1.0017     0.9983    0.9726     1.032
slope_prev       0.9789     1.0216    0.7642     1.254

Concordance= 0.97  (se = 0.015 )
Likelihood ratio test= 1.43  on 17 df,   p=1
Wald test            = 1.51  on 17 df,   p=1
Score (logrank) test = 1.52  on 17 df,   p=1

Linear models predicting mean differences

lm_meandiff_socdist <- lm(mean_diff_socdist ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_uk_cpt_socdist)
lm_meandiff_socdist %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_uk_cpt_socdist)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.59230 -0.47844 -0.01606  0.42775  1.91779 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -1.476e-15  6.443e-02   0.000  1.00000   
pers_o       1.733e-01  1.176e-01   1.474  0.14308   
pers_c      -1.721e-01  1.076e-01  -1.600  0.11221   
pers_e       3.030e-01  1.138e-01   2.662  0.00878 **
pers_a      -1.883e-01  9.469e-02  -1.988  0.04895 * 
pers_n      -1.439e-02  9.386e-02  -0.153  0.87842   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7375 on 125 degrees of freedom
Multiple R-squared:  0.4771,    Adjusted R-squared:  0.4562 
F-statistic: 22.81 on 5 and 125 DF,  p-value: 3.272e-16

lm_meandiff_socdist %>% confint(level=0.9)

                    5 %         95 %
(Intercept) -0.10677263  0.106772634
pers_o      -0.02157286  0.368159717
pers_c      -0.35034428  0.006189375
pers_e       0.11439435  0.491547403
pers_a      -0.34520521 -0.031373565
pers_n      -0.16992488  0.141149343

lm_meandiff_socdist_ctrl <- lm(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                            data = df_uk_cpt_socdist)
lm_meandiff_socdist_ctrl %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative + 
    onset_prev + slope_prev, data = df_uk_cpt_socdist)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.64573 -0.27449  0.03543  0.37169  1.53058 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.335186   0.291520   1.150 0.252680    
pers_o         0.038823   0.158645   0.245 0.807122    
pers_c        -0.134104   0.126954  -1.056 0.293095    
pers_e         0.074261   0.103326   0.719 0.473818    
pers_a         0.193685   0.100280   1.931 0.055955 .  
pers_n         0.121583   0.086393   1.407 0.162099    
airport_dist   0.055364   0.075845   0.730 0.466939    
males         -0.224919   0.089231  -2.521 0.013123 *  
popdens        0.579669   0.143703   4.034 0.000101 ***
manufacturing -0.232107   0.083830  -2.769 0.006587 ** 
tourism        0.093210   0.079389   1.174 0.242847    
health         0.012968   0.069975   0.185 0.853309    
academic      -0.020513   0.156637  -0.131 0.896043    
medinc         0.410681   0.103757   3.958 0.000133 ***
medage        -0.049986   0.137375  -0.364 0.716646    
conservative   0.162858   0.174256   0.935 0.352009    
onset_prev    -0.010575   0.009011  -1.174 0.243046    
slope_prev     0.038356   0.076465   0.502 0.616922    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6077 on 112 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.6797,    Adjusted R-squared:  0.6311 
F-statistic: 13.98 on 17 and 112 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl %>% confint(level=0.9)

                      5 %         95 %
(Intercept)   -0.14832058  0.818692212
pers_o        -0.22430100  0.301947876
pers_c        -0.34466665  0.076458268
pers_e        -0.09711227  0.245633955
pers_a         0.02736364  0.360006734
pers_n        -0.02170550  0.264871037
airport_dist  -0.07043049  0.181157779
males         -0.37291591 -0.076922591
popdens        0.34132736  0.818010783
manufacturing -0.37114475 -0.093069228
tourism       -0.03846205  0.224882948
health        -0.10309016  0.129026269
academic      -0.28030743  0.239281293
medinc         0.23859252  0.582770261
medage        -0.27783248  0.177860983
conservative  -0.12615763  0.451873581
onset_prev    -0.02551965  0.004369936
slope_prev    -0.08846702  0.165179585

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist <- cforest(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                               data = df_uk_cpt_socdist %>% drop_na(),
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp

       pers_o        pers_c        pers_e        pers_a        pers_n 
 0.0866398614  0.0139129422  0.0448894610  0.0243862299  0.0010570870 
 airport_dist         males       popdens manufacturing       tourism 
 0.0006984125  0.0063141561  0.2844528841  0.0896169809 -0.0019615256 
       health      academic        medinc        medage  conservative 
 0.0005118732  0.0287680803  0.1190340860  0.0127496646  0.0424411558 
   onset_prev    slope_prev 
 0.0024072003  0.0020732020 

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond

       pers_o        pers_c        pers_e        pers_a        pers_n 
 0.0827213621  0.0150431607  0.0415189305  0.0215636149  0.0009602196 
 airport_dist         males       popdens manufacturing       tourism 
 0.0008033941  0.0053098944  0.2567242510  0.0940480335 -0.0007021642 
       health      academic        medinc        medage  conservative 
 0.0009865419  0.0316531387  0.1111877822  0.0136577443  0.0394662118 
   onset_prev    slope_prev 
 0.0029093532  0.0004324219 

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist_2 <- coxph(Surv(cpt_day_socdist_2, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_uk_cpt_socdist)
cox_cpt_socdist_2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist_2, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_uk_cpt_socdist)

  n= 131, number of events= 131 

            coef exp(coef)  se(coef)      z Pr(>|z|)
pers_o -0.050902  0.950372  0.156351 -0.326    0.745
pers_c -0.060242  0.941537  0.149829 -0.402    0.688
pers_e  0.087182  1.091095  0.162827  0.535    0.592
pers_a  0.103047  1.108543  0.131924  0.781    0.435
pers_n  0.003533  1.003539  0.133928  0.026    0.979

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    0.9504     1.0522    0.6995     1.291
pers_c    0.9415     1.0621    0.7019     1.263
pers_e    1.0911     0.9165    0.7930     1.501
pers_a    1.1085     0.9021    0.8560     1.436
pers_n    1.0035     0.9965    0.7719     1.305

Concordance= 0.637  (se = 0.136 )
Likelihood ratio test= 0.92  on 5 df,   p=1
Wald test            = 0.93  on 5 df,   p=1
Score (logrank) test = 0.92  on 5 df,   p=1

# predict hazard from personality with controls
cox_cpt_socdist_ctrl_2 <- coxph(Surv(cpt_day_socdist_2, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                  data = df_uk_cpt_socdist)
cox_cpt_socdist_ctrl_2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist_2, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_dist + males + popdens + 
    manufacturing + tourism + health + academic + medinc + medage + 
    conservative + onset_prev + slope_prev, data = df_uk_cpt_socdist)

  n= 130, number of events= 130 
   (1 observation deleted due to missingness)

                    coef  exp(coef)   se(coef)      z Pr(>|z|)
pers_o        -0.1128691  0.8932676  0.2670282 -0.423    0.673
pers_c        -0.0342248  0.9663542  0.2083320 -0.164    0.870
pers_e         0.1041443  1.1097605  0.1768710  0.589    0.556
pers_a         0.0010809  1.0010814  0.1670260  0.006    0.995
pers_n        -0.0187305  0.9814438  0.1473178 -0.127    0.899
airport_dist  -0.0469830  0.9541036  0.1294395 -0.363    0.717
males          0.0373559  1.0380624  0.1579103  0.237    0.813
popdens       -0.0146463  0.9854604  0.2523392 -0.058    0.954
manufacturing  0.1499156  1.1617362  0.1336777  1.121    0.262
tourism        0.0253305  1.0256541  0.1308250  0.194    0.846
health         0.1119038  1.1184053  0.1163634  0.962    0.336
academic      -0.0233304  0.9769397  0.2648776 -0.088    0.930
medinc         0.0002381  1.0002381  0.1738722  0.001    0.999
medage         0.1011705  1.1064653  0.2282495  0.443    0.658
conservative  -0.1329131  0.8755412  0.3036219 -0.438    0.662
onset_prev     0.0010598  1.0010604  0.0153327  0.069    0.945
slope_prev     0.0194326  1.0196227  0.1327037  0.146    0.884

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           0.8933     1.1195    0.5293     1.508
pers_c           0.9664     1.0348    0.6424     1.454
pers_e           1.1098     0.9011    0.7847     1.570
pers_a           1.0011     0.9989    0.7216     1.389
pers_n           0.9814     1.0189    0.7353     1.310
airport_dist     0.9541     1.0481    0.7403     1.230
males            1.0381     0.9633    0.7617     1.415
popdens          0.9855     1.0148    0.6010     1.616
manufacturing    1.1617     0.8608    0.8940     1.510
tourism          1.0257     0.9750    0.7937     1.325
health           1.1184     0.8941    0.8903     1.405
academic         0.9769     1.0236    0.5813     1.642
medinc           1.0002     0.9998    0.7114     1.406
medage           1.1065     0.9038    0.7074     1.731
conservative     0.8755     1.1422    0.4829     1.588
onset_prev       1.0011     0.9989    0.9714     1.032
slope_prev       1.0196     0.9808    0.7861     1.323

Concordance= 0.788  (se = 0.127 )
Likelihood ratio test= 4.62  on 17 df,   p=1
Wald test            = 4.63  on 17 df,   p=1
Score (logrank) test = 4.64  on 17 df,   p=1

Linear models predicting mean differences

lm_meandiff_socdist_2 <- lm(mean_diff_socdist_2 ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_uk_cpt_socdist)
lm_meandiff_socdist_2 %>% summary()

Call:
lm(formula = mean_diff_socdist_2 ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_uk_cpt_socdist)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.22275 -0.46068 -0.02762  0.45220  1.96182 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.119e-15  5.933e-02   0.000  1.00000    
pers_o       2.500e-02  1.083e-01   0.231  0.81776    
pers_c      -1.594e-01  9.905e-02  -1.609  0.11018    
pers_e       4.230e-01  1.048e-01   4.037 9.39e-05 ***
pers_a      -2.424e-01  8.719e-02  -2.781  0.00627 ** 
pers_n      -1.862e-01  8.642e-02  -2.155  0.03309 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.679 on 125 degrees of freedom
Multiple R-squared:  0.5566,    Adjusted R-squared:  0.5389 
F-statistic: 31.39 on 5 and 125 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_2 %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.09831548  0.09831548
pers_o      -0.15442886  0.20443410
pers_c      -0.32350683  0.00478680
pers_e       0.24932419  0.59660402
pers_a      -0.38693082 -0.09795689
pers_n      -0.32944669 -0.04301177

lm_meandiff_socdist_ctrl_2 <- lm(mean_diff_socdist_2 ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                            data = df_uk_cpt_socdist)
lm_meandiff_socdist_ctrl_2 %>% summary()

Call:
lm(formula = mean_diff_socdist_2 ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative + 
    onset_prev + slope_prev, data = df_uk_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3539 -0.2995  0.0411  0.3338  1.1103 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.630291   0.234206   2.691 0.008211 ** 
pers_o        -0.023938   0.127455  -0.188 0.851359    
pers_c        -0.047208   0.101995  -0.463 0.644371    
pers_e         0.106542   0.083012   1.283 0.201980    
pers_a         0.094985   0.080565   1.179 0.240903    
pers_n        -0.036763   0.069408  -0.530 0.597387    
airport_dist  -0.095928   0.060934  -1.574 0.118239    
males         -0.248234   0.071688  -3.463 0.000758 ***
popdens        0.609775   0.115451   5.282 6.36e-07 ***
manufacturing -0.128026   0.067349  -1.901 0.059880 .  
tourism       -0.073726   0.063781  -1.156 0.250175    
health        -0.038582   0.056218  -0.686 0.493945    
academic       0.271148   0.125842   2.155 0.033330 *  
medinc         0.158529   0.083358   1.902 0.059770 .  
medage        -0.052320   0.110367  -0.474 0.636385    
conservative   0.220453   0.139997   1.575 0.118145    
onset_prev    -0.019925   0.007239  -2.752 0.006902 ** 
slope_prev     0.138325   0.061432   2.252 0.026293 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4882 on 112 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.7926,    Adjusted R-squared:  0.7612 
F-statistic: 25.18 on 17 and 112 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl_2 %>% confint(level=0.9)

                      5 %         95 %
(Intercept)    0.24184319  1.018739321
pers_o        -0.23533201  0.187455280
pers_c        -0.21637348  0.121957429
pers_e        -0.03113864  0.244222991
pers_a        -0.03863790  0.228606898
pers_n        -0.15188087  0.078354149
airport_dist  -0.19699032  0.005135194
males         -0.36713412 -0.129333676
popdens        0.41829143  0.801257933
manufacturing -0.23972830 -0.016322990
tourism       -0.17951095  0.032059895
health        -0.13182278  0.054659088
academic       0.06242939  0.479865921
medinc         0.02027300  0.296784706
medage        -0.23537132  0.130731893
conservative  -0.01174127  0.452647845
onset_prev    -0.03193196 -0.007918731
slope_prev     0.03643544  0.240214614

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist_2 <- cforest(mean_diff_socdist_2 ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                               data = df_uk_cpt_socdist %>% drop_na(),
                         controls = ctrls)

crf_meandiff_socdist_varimp_ <- varimp(crf_meandiff_socdist_2, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist_2, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp

       pers_o        pers_c        pers_e        pers_a        pers_n 
 0.0866398614  0.0139129422  0.0448894610  0.0243862299  0.0010570870 
 airport_dist         males       popdens manufacturing       tourism 
 0.0006984125  0.0063141561  0.2844528841  0.0896169809 -0.0019615256 
       health      academic        medinc        medage  conservative 
 0.0005118732  0.0287680803  0.1190340860  0.0127496646  0.0424411558 
   onset_prev    slope_prev 
 0.0024072003  0.0020732020 

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond

       pers_o        pers_c        pers_e        pers_a        pers_n 
 0.0105084218  0.0004282805  0.1288694275  0.0030573110  0.0016755109 
 airport_dist         males       popdens manufacturing       tourism 
 0.0579376000  0.0012207368  0.1855763980  0.0573411043  0.0032837033 
       health      academic        medinc        medage  conservative 
 0.0006151607  0.1741166028  0.0934911787  0.0053397776  0.0714661695 
   onset_prev    slope_prev 
 0.0058773245  0.0008061143 

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Export data

uk_list_results <- list(cox_onset_prev_laua, cox_onset_prev_laua_ctrl, 
     lm_slope_prev_laua, lm_slope_prev_laua_ctrl, 
     cox_cpt_socdist, cox_cpt_socdist_ctrl,
     lm_meandiff_socdist, lm_meandiff_socdist_ctrl,
     cox_cpt_socdist_2, cox_cpt_socdist_ctrl_2,
     lm_meandiff_socdist_2, lm_meandiff_socdist_ctrl_2)

results_names <- list('cox_onset_prev', 'cox_onset_prev_ctrl', 
     'lm_slope_prev', 'lm_slope_prev_ctrl', 
     'cox_cpt_socdist', 'cox_cpt_socdist_ctrl', 
     'lm_meandiff_socdist', 'lm_meandiff_socdist_ctrl',
     'cox_cpt_socdist_2', 'cox_cpt_socdist_ctrl_2',
     'lm_meandiff_socdist_2', 'lm_meandiff_socdist_ctrl_2')

names(uk_list_results) <- results_names

save(uk_list_results, file="uk_list_results.RData")

write_csv(df_uk_slope_prev, '/Users/hp2500/Google Drive/STUDY/Columbia/Research/Corona/Delivery/df_uk_slope_prev.csv')
write_csv(df_uk_cpt_socdist, '/Users/hp2500/Google Drive/STUDY/Columbia/Research/Corona/Delivery/df_uk_cpt_socdist.csv')